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I.  INTRODUCTION 


In  situations  where  the  full  two-fluid  plasma  system(23)  is  difficult  to  solve,  near  the 
MHD  limit  for  example  where  the  Hall  parameter  is  small,  a  simpler  approach  is  to  solve 
the  two-temperature  Hall  MHD  model  with  diamagnetic  drift.  This  model  eliminates  the 
displacement  current  as  well  as  charge  separation  and  the  electron  cyclotron  frequency  while 
still  capturing  many  two-fluid  effects.  Hall  MHD  has  been  studied  by  many  authors  and 
several  algorithms  have  been  developed  (3;  11;  13;  14).  The  discontinuous  Galerkin  technique 
has  been  applied  to  many  plasma  systems  (16;  17;  20;  24)  and  the  technique  is  discussed  in 
great  detail  in  general  terms  in  the  textbook  by  Hesthaven  and  Warburton  (12).  James  Wiley 
of  University  of  Texas  at  Austin  also  investigated  DG  method  for  the  Hall  MHD  system 
though  that  work  remains  un-published.  This  work  is  unique  in  that  the  discontinuous 
Galerkin  method  is  used  and  the  algorithm  is  developed  from  the  Hall  MHD  equations  in 
gas  dynamic  conservative  form.  It  is  believed  that  this  form  is  appropriate  because  the 
magnetic  held  and  ion  fluid  decouple  in  regimes  where  the  Hall  term  is  important  and 
therefore  shocks  appear  that  obey  the  ion  jump  conditions  only.  Currently  these  results  are 
restricted  to  second  order  spatial  accuracy  and  3rd  or  4th  order  temporal  accuracy. 

II.  SYSTEM  OF  EQUATIONS 

The  model  is  derived  by  assuming  that  the  speed  of  light  is  large  in  comparison  to  other 
speeds  in  the  system.  Quasi-neutrality  is  assumed,  so  that  only  one  continuity  equation  is 
required, 

^  +  V|p,C/,]=0,  (1) 

Along  with  a  momentum  equation  for  the  ions, 

^l  +  s/.\PiuiUl  +  Pi\  =  nlqlE  +  Jt.xB1  (2) 

an  ion  energy  equation, 

Qe . 

-^  +  \/-[Ul(el  +  Pl)\  =  Jl-E,  (3) 
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where  ex  is  the  ion  kinetic  energy  +  internal  energy,  and  an  electron  energy  equation, 


Qe 

+  V  ■[Ue(ee  +  Pe)\  =  Je-E. 


(4) 


where  ee  is  the  electron  kinetic  energy  +  the  electron  internal  energy  -  however,  electron 
kinetic  energy  is  zero  because  the  electron  mass  is  set  to  zero.  The  species  currents  are  Ja 
where  a  —  e,i  for  the  electrons  and  ions  respectively.  The  ideal  gas  equation  of  state  is 
used,  therefore  P*  =  rq/cP;  and  Pe  =  nekTe.  The  completeness  of  the  system  requires  a 
momentum  equation  for  the  electron  fluid.  Since  the  electrons  are  assumed  massless  the 
convective  derivative  is  set  to  zero  and  we  are  left  with  pressure  and  Lorentz  forces.  The 
reduced  electron  momentum  equation  or  Ohm’s  law,  is  given  by 


T  —  neqeE  —  —  Je  x  B  +  V  Pe 


(5) 


The  magnetic  held  equation  is  also  included 


dB  _  ^ 

_  +  Yx£  =  0 


(6) 


In  this  paper  divergence  cleaning  is  accomplished  using  the  hyperbolic  divergence  cleaning 
method  as  suggested  in  (9).  Using  this  formulation  an  additional  equation  is  added  and  the 
magnetic  held  equation  is  modified,  thus 


dB 

— -  +  VxP  +  V0  =  0 
at 


(7) 


and 


^  +  r2v  ■  b  = -«. 


(8) 


The  modified  magnetic  held  equation  adds  the  divergence  wave  to  the  system  which  prop¬ 
agates  divergence  errors  out  of  the  domain  at  speed  T  while  simultaneously  damping  them. 
A  linear  wave  analysis  shows  that  this  divergence  wave  does  not  couple  to  the  plasma  and 
therefore  does  not  modify  any  other  plasma  wave. 
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III.  DISPERSION  RELATION 


For  complicated  systems  the  wave  systems  may  contain  dispersive  waves  where  the  fre¬ 
quency  lv  depends  on  wave  number  k  in  a  non-linear  fashion.  Hall  MHD  in  particular, 
admits  solutions  where  u  is  quadratic  in  k.  A  careful  analysis  of  dispersion  relations  is  often 
necessary  in  order  to  properly  determine  time  step  restrictions.  The  following  can  be  said 
about  Hall  MHD, 

The  Hall  MHD  model  produces  the  whistler  wave  in  the  case  where  the  wave  travels 
parallel  to  the  magnetic  held 

B0  2  B(i  \J 4  n0  rrii  /j0  q2  +  k2  rnf 

w= - k  H - - - -  k  .  (9) 

2  qn0n0  2qmin0fi0 

Note  that  this  can  be  written 


w  =  ^uci  d2  k'2  +  J V\k 2  +  i  u2ci  d\  kA , 


where  d,  = 


JL.  =  CJ  M  i 


is  the  ion  inertial  length  or  ion  skin  depth,  Va  =  ,  Bo 

n0qz  &  r  i  A  y/nom.i  /xq 


is  the 


Alfven  wave  speed  and  Vw  =  c jrjd?  =  B(<  is  the  whistler  coefficient.  The  whistler  wave 

^  w  c  1  i  q  no  /j,o 

becomes, 

u  =  uci  d 2  k2  (11) 

for  sufficiently  large  k.  The  importance  of  the  Hall  term  (and  hence  the  whistler  wave)  can 
be  determined  by  the  ratio  of  the  contributions  of  the  quadratic  and  linear  terms  in  the 
whistler  wave  dispersion  relation,  thus,  the  whistler  wave  is  important  when, 


’-k  =  d,k  A  1 


This  condition  implies  that  the  Hall  term  is  important  for  frequencies  where  u  ^  u)Ci  and 
for  wavelengths  L  =  ^  where  L  ;$  d{. 

For  small  wave  numbers  the  whistler  wave  becomes  the  Alfven  wave  and  u  varies  linearly 
with  k.  In  the  limit  of  long  wave  lengths  u  varies  quadratically  with  k,  thus  we  have  tradi¬ 
tional  convection  combined  with  dispersive  convection.  This  sort  of  consideration  suggests 
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the  following  system  as  a  simple  model  for  the  ideal  Hall  MHD  system 


d  u 
~dt 


r  du  •  r 

m  TT - 102 

OX 


d2u 
dx2  ' 


(13) 


If  the  non-dispersive  advective  term  is  ignored  we  simply  have, 


du  d  2u 

~di  =  ~~ 162  d 


(14) 


which  is  equivalent  to  the  Schrodinger  equation.  Fortunately,  numerical  solutions  to  the 
Schrodinger  equation  have  been  studied  for  decades,  see  for  example  (4),  and  basic  stable 
schemes  can  be  borrowed  from  previous  work. 

In  the  case  where  the  wave  travels  perpendicular  to  the  magnetic  field  the  dispersion 
relation  produces  the  magnetosonic  wave 

Bl  +  Pq7/x0 

TTli  TIq  /i 0 

which  is  non-dispersive  and  linear  in  k. 


A.  Time  Step  Restrictions 


Using  an  explicit  code  the  following  time  step  restrictions  are  suggested.  First  of  all  the 
magnetosonic  wave  Vm  should  be  resolved,  thus 


At  < 


Ax 


\u\  +  vm 


(16) 


From  the  whistler  wave  dispersion  relation  the  maximum  time  step  should  be 


At  < 


Ax2 

"kT 


In  general  then, 


At  <  min 


Ax  Ax2 


(17) 


(18) 


[\U\  +  VJ  Vw  J 

Drift  waves  also  exist  in  the  Hall  MHD  system  and  the  actual  time  step  may  be  more 
restrictive.  Huba  discusses  this  in  more  detail  in  (14).  For  many  problems  the  Whistler 
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wave  dominates  and  taking  a  time  step  based  on  its  value  is  frequently  adequate. 


IV.  NUMERICAL  APPROACH 

The  numerical  approach  used  is  the  discontinuous  Galerkin  method.  The  discontinuous 
Galerkin  method  can  be  considered  an  extension  of  finite  volume  methods  or  a  finite  element 
method  where  the  solution  can  be  discontinuous  at  nodes.  The  technique  is  described  for 
hyperbolic  systems  in  (5-8).  The  method  is  explicit  and  does  not  require  the  solution  of  a 
global  mass  matrix  at  each  time  step. 

Start  by  defining  n  auxiliary  variables.  In  the  case  of  Hall  MHD  with  diamagnetic  drift, 
n  is  2.  These  auxiliary  variables  are  written  generically  as, 


V\  =  V  ■  Wi  (q,  Vi)  +  (q,  Vj,r)  . 


(19) 


't>2  =  V  •  w2  (q,  t>2,  t>i)  +  (q,  t>2,  Vi,  r)  . 


(20) 


and  finally, 


v,t  =  V  ■  wn  (q,vn- +'tpn  (q, ...,«i,r)  . 


(21) 


The  q  are  “conserved”  variables,  the  v  are  auxiliary  variables,  the  w  are  flux  functions  and 
the  "0  are  terms  that  only  depend  on  algebraic  combinations  of  q,  v,  and  r  where  r  is  a 
vector  position.  Using  the  above  formulation  the  Hall  MHD  system  can  be  described.  The 
auxiliary  variables  for  the  Hall  MHD  system  are  defined  as, 


V\  —  J  —  — V  x  B 
/U 

and 

v2  =  'b  =  ne  qe  E  =  -  ( J  -  Ji)  x  B  +  V  Pe  . 
In  addition,  take  the  balance  law, 


dq 

~dt 


+  V  ■  f  (q,vn,,..,v2,v  i)  =  i>(q,  li2,  tii,  r) 


(22) 


(23) 


(24) 


where  /  is  a  flux  function. 
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The  balance  law  represents  each  of  the  equations  with  time  derivatives  in  them  in  the 
Hall  MHD  system.  The  approach  is  to  solve  for  the  auxiliary  variables  one  after  another, 
using  the  solution  un_i  to  solve  for  vni  then  substitute  these  solution  into  the  balance  law. 

Using  quadrilateral  elements,  Legendre  polynomials  are  used  as  basis  function  and  each 
of  the  variables  vn  and  qn  are  represented  as  a  linear  combination  of  Legendre  polynomials. 
Therefore  if  the  basis  set  is  hr,  where  hr  is  the  rth  Legendre  polynomial,  then  the  variables 
are  represented  as, 

p- 1 

vn  ^  ^  oinrhr  ,  (25) 

r=0 

v- 1 

Qn  ^  ^  finrhr  •>  (26) 

r= 0 

where  p  is  the  number  of  basis  functions  in  the  expansion.  In  1-D,  p  corresponds  to  the 
order  of  accuracy  of  the  scheme.  In  2-D  for  a  strictly  nth  order  scheme  p  =  \  (n  +  1)  n 
where  2D  Legendre  polynomials  are  used  as  the  bases.  In  equation  (25)  anr  represents  the 
rth  basis  coefficient  of  the  nth  auxiliary  variable  and  in  equation  (26)  (3nr  represents  the  rth 
coefficient  of  the  nth  “conserved”  variable.  For  a  system  of  one  auxiliary  variable  and  one 
conserved  variable  you  can  eliminate  the  subscript  n  to  see  that  each  auxiliary  variable  and 
each  conserved  variable  is  simply  a  sum  of  Legendre  polynomials.  These  same  basis  functions 
are  used  to  expand  the  system  of  equations.  First  the  auxiliary  variables  are  considered. 
Multiply  the  equation  (21)  by  the  basis  function  set  hr  and  integrate  over  the  volume  to  get, 


vn  hr  dfl  = 


V  •  wn  (, q ,  vn—  i , ...,  Ui)  hr  dQ  +  /  i\)n  (q,  un_i, ...,  vu  r)  hr  dtt  (27) 


'n 


Integration  by  parts  is  applied  to  the  derivative  to  yield, 


vn  hr  dil  —  /  [wn  (q,  un_i, ...,  tq)  hr\  ■  n  ddVt— 

Jon 

[  wn(q,vn-i,...,v1)-\7hrdQ+  [  ipn  (q,  un_i, ...,  vu  r)  hr  dfl  =  Tnr  (28) 


Substituting  the  definition  of  vn  (25)  into  the  first  term  of  equation  (28),  the  following  result 
is  obtained 


(29) 
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so  that 


(30) 


A 


gr 


hg  hr  dit . 


Jn 

On  regular  grid,  Agr  is  diagonal  since  the  basis  functions  are  orthogonal.  Now  the  coefficients 
ang  can  be  calculated  based  on  a  numerical  solution  to  the  right  hand  side  of  (28),  so  that 


O^ng  — 


Tn  r  [A 


1-1 


gr  J 


(31) 


Agr  is  a  square  matrix  of  size  p  x  p,  ang  is  a  matrix  of  size  m  x  p  where  m  is  the  number  of 
auxiliary  variables  in  the  system  and  Tnr  is  a  m  x  p  matrix.  The  right  hand  side  is  solved 
by  numerically  integrating  the  integrals  at  each  time  step.  Agr  is  evaluated  at  the  beginning 
of  the  simulation.  The  only  derivatives  that  remain  are  on  the  basis  function  hr  and  those 
can  be  calculated  either  analytically  in  the  case  of  cartesian  coordinates  or  numerically  for 
general  geometries.  The  terms  that  remain  wn  (q,  u„_i, ...,  iq)  and  V'n  (?,  vn-i,  •••,  ffi,  r)  are 
algebraic  function  of  q  and  vn  which  are  known  (as  a  result  of  the  basis  function  expansion 
of  equations  (25)  and  (26))  at  all  points  in  space.  Gaussian  quadrature  is  performed  to 
approximate  the  integrals.  For  a  second  order  method  in  ID  a  2  point  gaussian  quadrature 
is  performed.  In  general  in  n  dimensions  a  pn  point  gaussian  quadrature  is  performed  for 
an  order  p  accurate  solution. 

The  balance  laws  are  solved  in  a  virtually  identical  manner  except  a  time  derivative 
must  also  be  evaluated,  for  completeness,  start  with  the  balance  law  multiplied  by  the  basis 
function  set  hr . 


dqn 

dt 


hrdQ+  /  V-f(q,v  i,v2,...,vn)  hrdU=  /  (q,  vn, ...,  v2,  v1,  r)  hrdtt. 


(32) 


move  the  second  term  to  the  right  hand  side  and  expand  the  derivative  using  integration  by 
parts 


d  q. 


—  hr  dQ  = 
dt 


'an 


[f  (q,vn,...,v2,v i)  hr]-nddQ+ 


f  (q,vn,...,v2,v i)  •  \/hrdQ+  /  'ip(q,vn,...,v2,v1,r)  hrdtt  =Ynr  (33) 


The  definition  of  q  is  substituted  into  the  first  term  so  that  the  first  term  becomes 


This  gives  us  equations  for  the  time  derivative  of  /3ng 

^  =  Ynr\Alr}-\  (35) 

Once  the  right  hand  side  of  equation  (35)  is  calculated  numerically,  the  (3ng  can  be  calculated 
using  a  standard  second  order  Runge-Kutta  technique.  Note  that  (3ng  is  a  mxp  matrix  and 
so  is  Ynr.  Notice  that  the  right  hand  side  of  equation  (33)  is  identical  in  form  to  (28).  The 
right  hand  side  of  (33)  is  therefore  calculated  in  an  identical  manner. 

A.  Numerical  Fluxes 

One  important  issue  in  applying  discontinuous  Galerkin  methods  is  the  evaluation  of  the 
numerical  flux.  At  boundary  interfaces  the  flux  is  undefined  because  the  solution  is  generally 
discontinuous  at  that  point.  A  simple  flux,  known  as  the  local  Lax  Friedrichs  flux  can  be 
used  on  the  conserved  variables,  in  this  case 

F,+ 1/2  =  \  (Ft  +  -  ||A|  -  q?)  (36) 

in  general,  for  auxiliary  variables  this  same  flux  can  be  used,  but  it  needs  to  be  modified  by 
the  factor  At,  namely 

F<+ !/2  =  \  (Ft  +  F-+1)  -  iAt|A|  (q-+1  -  qf)  (37) 

Choice  of  fluxes  for  diffusive  systems  are  discussed  in  (15)  and  flux  functions  for  the 
Schrodinger  equation  are  discussed  in  (19).  Both  papers  indicate  that  centered  fluxes  will  be 
stable,  the  added  dissipation  may  help  in  avoiding  negative  pressures  in  certain  situations. 
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V.  RESULTS 


A.  Brio  and  Wu  Shock 

The  Brio  and  Wu  shock  problem  was  recently  extended  to  systems  with  finite  ion  inertial 
length  effects  in  the  following  works  (10;  16;  23)  in  the  framework  of  the  two-fluid  plasma 
system.  These  same  simulations  can  be  run  using  Hall  MHD  with  diamagnetic  drift.  The 
results  that  follow  are  Hall  MHD  solutions  to  these  two-fluid  shock  problems. 

Parameters  used  are  qt  =  10,  qe  =  —10,  and  e0  =  1,  fj,0  =  1,  c  =  1,  ye  =  7 \  —  |, 
rrii  =  1.  The  initial  conditions  are  given  by  Pe  =  Pi  =  0.5  x  10”4,  ne  =  nt  =  1.0  and 
B  =  (7.5  x  10~3, 1.0  x  10~2,  0)  on  the  left  half  of  the  domain.  On  the  right  half  of  the 
domain  the  Pe  =  Pt  =  0.5  x  10~5,  ne  =  rii  =  0.125  and  B  =  (7.5  x  10~3,  -1.0  x  10~2,0) 
while  Ue  =  Ui  =  E  =  0  everywhere. 

Figures  1,2  and  3  show  the  Brio  and  Wu  shock  solution  with  rgi  =  1/10,  rgi  =  1/100  and 
rgi  =  1/1000  respectively  where  rgi  is  the  ion  gyroradius  given  the  initial  conditions  on  the 
left  hand  side  of  the  shock.  The  ion  gyroradius  is  calculated  as  vthi/uCi  where  vthi  is  the  ion 
thermal  velocity  and  u oCi  =  Solutions  are  quite  similar  to  those  calculated  previously 

using  the  two-fluid  equations (10). 

The  two-fluid  results  are  virtually  identical  to  the  Hall  MHD  results  so  it’s  tempting  to 
suggest  that  Hall  MHD  should  be  used  in  the  vast  majority  of  cases  because  the  physics  is 
simpler.  Difference  only  begin  to  appear  at  the  Debye  length.  In  regimes  where  the  speed  of 
light  is  high  compared  to  fluid  speeds  the  speed  of  light  can  be  artificially  reduced  in  the  full 
two-fluid  system.  The  reduced  speed  of  light  and  electron  inertia  act  to  provide  numerical 
limits  for  dispersive  waves  thus  reducing  the  stiffness  of  the  system  compared  to  Hall  MHD. 
The  displacement  current  and  electron  inertia  can  be  thought  of  as  numerical  parameters 
that  would  reduce  the  stiffness  of  the  Hall  MHD  system  just  as  artificial  viscosity  is  often 
used  to  stabilize  numerical  solutions  to  purely  hyperbolic  otherwise  non  viscous  system.  The 
full  two-fluid  system  can  thus  be  thought  of  as  the  Hall  MHD  system  with  numerical  terms 
that  limit  the  stiffness  of  Hall  MHD.  When  actually  running  the  simulations  in  the  regime 
where  rgi  «  1/10  the  domain  or  larger,  the  full  two-fluid  simulations  complete  considerably 
faster  than  the  Hall  MHD  simulations.  At  the  other  extreme,  where  rgi  ~  1/1000  the  Hall 
MHD  simulations  complete  considerably  faster  than  the  two-fluid  simulations  provided  that 
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Mass  density  with  characteristic  ion  Larmor  radius  of  1/10 


x 


FIG.  1  Hall  MHD  solution  to  the  Brio  and  Wu  shock  problem  with  characteristic  ion  Larmor  radius 
of  1/10.  The  Hall  MHD  solution  agrees  quite  well  with  published  two-fluid  solutions.  In  this  case 
the  obvious  differences  are  Debye  length  oscillations  which  occur  at  the  shock  (and  elsewhere)  in 
the  two-fluid  solution. 


Mass  density  with  characterstic  ion  inertial  length  of  1/100 


FIG.  2  Hall  MHD  solution  to  the  Brio  and  Wu  shock  problem  with  characteristic  ion  inertial 
length  1/100.  The  Hall  MHD  solution  agrees  quite  well  with  published  two-fluid  solutions.  Major 
differences  include  the  large  oscillations  on  the  left  of  the  shock.  These  oscillation  also  appear  in 
the  two-fluid  solutions  when  sufficiently  small  time  step  are  taken. 
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Mass  density  with  characterstic  ion  inertial  length  of  1/1000 


FIG.  3  Hall  MHD  solution  to  the  Brio  and  Wu  shock  problem  with  characteristic  ion  inertial  length 
1/1000.  The  solution  agrees  remarkably  well  with  published  two-fluid  solutions.  In  this  regime 
the  two-fluid  solution  is  particularly  difficult  to  solve  due  to  the  stiffness  of  the  system.  The  Hall 
MHD  solution  on  the  other  hand  is  quite  simple  because  the  solution  is  a  small  perturbation  from 
the  MHD  solution. 

rgi  is  not  well  resolved.  Moving  towards  the  MHD  limit  where  the  whistler  wave  is  not 
well  resolved,  solving  the  Hall  system  is  still  computationally  efficient,  however  as  resolution 
increases,  the  whistler  wave  speed  increases  and  a  technique  that  limits  the  whistler  wave 
speed  is  desirable.  The  two-fluid  model  becomes  computationally  inefficient  near  the  MHD 
limit  because  cyclotron  and  plasma  frequency  time  scales  must  be  resolved  along  with  Debye 
length  spatial  scales  so  that  much  larger  computational  grids  are  necessary  than  are  required 
for  Hall  MHD  where  these  same  time  and  spatial  scales  do  not  need  to  be  resolved. 

Relative  computational  times  can  be  estimated  by  looking  at  the  waves  involved.  The 
full  two-fluid  model  can  be  solved  in  roughly  the  same  manner  as  the  Hall  MHD  system, 
so  one  step  of  the  Hall  MHD  algorithm  takes  about  as  long  as  one  step  of  a  full  two-fluid 
algorithm.  Therefore,  simulation  time  depends  on  the  necessary  spatial  resolution  and  time 
step  restriction.  In  full  two-fluid  At  =  A x/c  where  c  is  the  speed  of  light.  For  Hall  MHD 
the  fastest  wave  is  given  approximately  by  ujcldf  /Ax  so  that  At  =  Ax\ .  Computational 
effort  required  for  each  system  is  given  roughly  by  A^At  so  that  the  ratio  of  computational 
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effort  (time)  is 


c  Ax 


(38) 


Etf 

E}Jall  ^  c  i  d'j 

where  Etf  is  the  computation  time  required  to  solve  the  full  two-fluid  system  at  the  given 
resolution  and  Euaii  is  the  computational  time  required  to  solve  the  Hall  MHD  system  at  the 
same  resolution.  It  is  clear  that  as  resolution  increases,  i.e.  Ax  decreases,  the  full  two-fluid 
system  requires  less  time  to  solve  relative  to  the  Hall  MHD  system. 

Ultimately  one  could  argue  that  implicit  schemes  are  of  greater  importance  and  time 
constraints  are  less  dependent  on  grid  resolution  in  that  case,  however  the  above  also  suggest 
that  at  high  resolution  the  full  two-fluid  system  is  less  stiff  than  the  Hall  MHD  system  and 
faster  convergence  might  be  observed  using  implicit  methods  on  the  full  two-fluid  system 
than  on  the  Hall  MHD  system. 
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B.  Magnetic  Reconnection 


Magnetic  reconnection  is  an  important  test  of  Hall  MHD  codes  and  has  become  a  stan¬ 
dard  benchmark  of  Hall  MHD  and  Two-Fluid  codes.  The  following  results  illustrate  Hall 
reconnection  with  results  compared  to  solutions  computed  using  other  methods.  A  good  ref¬ 
erence  for  collisionless  reconnection  can  be  found  in  (2),  and  in  the  GEM  challenge  papers 
(21;  22). 

The  GEM  challenge  magnetic  reconnection  problem  is  non-dimensionalized  as  in  (1) 

where  lengths  are  normalized  by  the  ion  inertial  length  d  =  c/wpi  =  c  j  2  time  is  non- 

dimensionalized  by  the  ion-cyclotron  time  where  B0  is  the  magnetic  held  at  infinity.  The 

1 

(  B ^  \  2 

velocities  are  normalized  by  the  Alfven  velocity  Vn  =  - u —  .  Finally  current  density 

•7  J  u  \  no  rrii  no  J  J  J 

is  non-dimensionalized  by  Jo  =  B°Wpi  and  E  by  E0  =  VaB0.  The  domain  is  (— 6.4  d,  6.4  d) 
and  the  simulation  is  run  out  to  40/wc;.  Conducting  walls  are  used  on  the  y  boundaries 
and  periodic  boundaries  are  used  on  the  x  boundaries.  A  =  0.5  d  and  the  specihc  heat  ratio 
7  =  |.  The  initial  number  densities  are  given  by, 

n  =  n0  Q  +  sech2  .  (39) 

The  electron  and  ion  temperatures  differ  slightly  from  the  GEM  challenge  problem  in  that 
the  electron  and  ion  temperatures  are  equal,  this  gives  the  following  electron  pressure,  Pei 


Pe 


- Bn  — 

12/7,0  u  n0 


and  ion  pressure  Pi 


R  =  —Bl 


12/^0  no 

The  electron  and  ion  pressure  balance  the  magnetic  held  which  is  given  by 


d  d  +.  i  (y\  L  B° 71  ( 2itx\  ■  ( ny 

B*  =  Bo  tanh  ( A )  +  To  X 008  (— ) 3111  (y; 


B 


y 


Bo 
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(40) 


(41) 


(42) 

(43) 
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GEM  reconnection  using  Two-Fluid  and  Hall  MHD 


FIG.  4  Plots  of  reconnected  flux  calculated  using  the  full  two-fluid  system  and  the  Hall  MHD 
system  discussed  in  this  paper. 


The  magnetic  field  is  in  equilibrium  with  with  the  current  J, 

J,e  =  ^sech20  .  (44) 

However,  for  the  Hall  MHD  system,  J  is  calculated  numerically  from  B  so  it  is  not  necessary 
to  set  it  as  an  initial  condition. 

Figure  4  shows  the  reconnected  flux  vs.  time  for  the  Hall  MHD  solution  compared  to  the 
two-fluid  solution  calculated  in  (16).  The  reconnection  rates  agree  fairly  well  as  expected, 
both  solutions  produce  fast  reconnection.  Figure  5  shows  plots  of  total  current  using  the 
Hall  MHD  code  at  time  t  =  28/u;Cj.  The  Hall  MHD  solution  shows  much  less  structure  than 
the  two-fluid  solution;  this  is  likely  a  result  of  the  neglect  of  electron  inertia,  the  electron 
convective  derivative  in  particular,  which  allows  shear  instabilities  in  the  electron  fluid  in 
the  two-fluid  case. 
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FIG.  5  A  Low  resolution  simulation  of  Hall  reconnection.  Growth  rate  is  somewhat  lower  than 
expected  for  a  high  resolution  simulation. 


C.  Axisymmetric  Z-Pinch 

The  Z-pinch  is  a  good  test  of  the  axisymmetric  code.  In  this  case  a  new  benchmark  is 
presented  based  on  the  Hall  solution  to  the  axisymmetric  Z-pinch  originally  presented  in 
(18)  which  used  the  full  two-fluid  model.  The  Hall  solution  is  important  because  now  the 
two  techniques  can  be  compared  and  differences  in  the  two  plasma  models  identified.  Both 
Hall  MHD  and  the  full  two-fluid  model  show  the  development  of  a  fast  instability  on  the 
scale  of  the  ion  inertial  length  dj.  In  the  Hall  MHD  solution  shown  in  figure  6  the  structure 
is  a  finer,  this  may  be  explained  by  the  fact  that  in  Hall  MHD  the  electron  mass  is  zero  so 
electron  inertia  effects  are  ignored. 
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FIG.  6  Mass  density  of  a  Z-pinch  where  the  axis  are  in  units  of  ion  Larmor  radii.  The  solution 
shows  the  growth  of  a  fast  instability  much  like  that  in  Figure  8  of  (18)  where  the  same  initial 
conditions  are  used. 
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VI.  CONCLUSION 


A  Hall  MHD  code  including  diamagnetic  drift  has  been  developed  using  the  discontinuous 
Galerkin  method.  The  technique  is  explicit  and  involves  a  straight  forward  application  of 
the  DG  method.  In  principal  this  same  algorithm  can  be  extended  to  3rd  and  higher  order 
accuracy  and  can  be  applied  to  non-Cartesian  grids.  The  simulations  of  Brio  and  Wu 
shock  show  good  agreement  with  those  performed  using  the  full  two-fluid  model,  however 
it  has  been  observed  that  the  full  two-fluid  model  is  actually  achieves  solutions  faster  than 
Hall  MHD  in  regimes  where  the  Hall  parameter  is  large  compared  to  the  domain  size. 
Simulations  of  magnetic  reconnection  agree  with  prior  published  results  as  do  simulations 
of  fast  instabilities  in  Z-pinches. 
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